Micron Document




Formule BBP
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
top
La formule BBP (ou formule de Bailey-Borwein-Plouffe) permet de calculer le n-ième chiffre après la virgule du nombre π en base 2 (ou 16) sans avoir à en calculer les précédents, et en utilisant très peu de mémoire et de temps. Elle a été obtenue le 19 septembre 1995 par Simon Plouffe en collaboration avec David H. Bailey et Peter Borweincite-ref-1[1].

Contents


──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

La formule

Dans sa forme originelle, la formule BBP est donnée par

π π = ∑ ∑ k = 0 ∞ ∞ 1 16 k ( 4 8 k + 1 − − 2 8 k + 4 − − 1 8 k + 5 − − 1 8 k + 6 ) . {\displaystyle \pi =\sum _{k=0}^{\infty }{\frac {1}{16^{k}}}\left({\frac {4}{8k+1}}-{\frac {2}{8k+4}}-{\frac {1}{8k+5}}-{\frac {1}{8k+6}}\right).}

Démonstration

On note S n = ∑ ∑ k = 0 ∞ ∞ 1 16 k ( 8 k + n ) {\displaystyle S_{n}=\sum _{k=0}^{\infty }{\frac {1}{16^{k}(8k+n)}}} et on démontre ici la formule de Plouffe généralisée :

( 0 ) ∀ ∀ r ∈ ∈ C π π = ( 4 + 8 r ) S 1 − − 8 r S 2 − − 4 r S 3 − − ( 2 + 8 r ) S 4 − − ( 1 + 2 r ) S 5 − − ( 1 + 2 r ) S 6 + r S 7 {\displaystyle (0)\quad \forall r\in \mathbb {C} \qquad \pi \ =\ (4+8r)S_{1}-8rS_{2}-4rS_{3}-(2+8r)S_{4}-(1+2r)S_{5}-(1+2r)S_{6}+rS_{7}} .
(le cas r = 0 est sa formule originale ; le cas r = –1/4 est, sous une forme plus détaillée, celle d'Adamchick et Wagon).

On pose α α = 1 − − i = 2 e − − i π π / 4 {\displaystyle \alpha =1-\mathrm {i} ={\sqrt {2}}\ \mathrm {e} ^{-\mathrm {i} \pi /4}} et on va calculer de deux façons l'intégrale suivante :

I

=



0

1

d

y

α

α



y

.

{\displaystyle I=\int _{0}^{1}{\frac {\mathrm {d} y}{\alpha -y}}.}

Elle est d'une part reliée aux Sn par :

I = 1 α α ∫ ∫ 0 1 d y 1 − − y / α α = 1 α α ∫ ∫ 0 1 ∑ ∑ m = 0 ∞ ∞ y m α α m d y = ∑ ∑ m = 0 ∞ ∞ e i ( m + 1 ) π π / 4 ( m + 1 ) 2 m + 1 = 1 + i 2 S 1 + i 2 S 2 + − − 1 + i 4 S 3 − − 1 4 S 4 − − 1 + i 8 S 5 − − i 8 S 6 + 1 − − i 16 S 7 + 1 16 S 8 , {\displaystyle {\begin{aligned}I&={\frac {1}{\alpha }}\int _{0}^{1}{\frac {\mathrm {d} y}{1-y/\alpha }}={\frac {1}{\alpha }}\int _{0}^{1}\sum _{m=0}^{\infty }{\frac {y^{m}}{\alpha ^{m}}}\mathrm {d} y=\sum _{m=0}^{\infty }{\frac {\mathrm {e} ^{\mathrm {i} (m+1)\pi /4}}{(m+1){\sqrt {2}}^{m+1}}}\\&={\frac {1+\mathrm {i} }{2}}S_{1}+{\frac {\mathrm {i} }{2}}S_{2}+{\frac {-1+\mathrm {i} }{4}}S_{3}-{\frac {1}{4}}S_{4}-{\frac {1+\mathrm {i} }{8}}S_{5}-{\frac {\mathrm {i} }{8}}S_{6}+{\frac {1-\mathrm {i} }{16}}S_{7}+{\frac {1}{16}}S_{8},\end{aligned}}}

et d'autre part calculable par des méthodes élémentaires (en calculant séparément sa partie réelle et sa partie imaginaire), ou de façon plus synthétique via le logarithme complexe :

I = − − [ ln ⁡ ⁡ ( α α − − y ) ] 0 1 = ln ⁡ ⁡ ( α α α α − − 1 ) = ln ⁡ ⁡ ( 1 + i ) = ln ⁡ ⁡ ( 2 e i π π / 4 ) = ln ⁡ ⁡ 2 2 + i π π 4 . {\displaystyle {\begin{aligned}I&=-\left[\ln(\alpha -y)\right]_{0}^{1}=\ln \left({\frac {\alpha }{\alpha -1}}\right)=\ln(1+\mathrm {i} )=\ln({\sqrt {2}}\ \mathrm {e} ^{\mathrm {i} \pi /4})={\frac {\ln 2}{2}}+\mathrm {i} {\frac {\pi }{4}}.\end{aligned}}}

L'égalité entre ces deux expressions de I équivaut à :

( 1 ) π π = 4 I m ( I ) = 2 S 1 + 2 S 2 + S 3 − − 1 2 S 5 − − 1 2 S 6 − − 1 4 S 7 , {\displaystyle (1)\qquad \pi =4\ \mathrm {Im} (I)=2S_{1}+2S_{2}+S_{3}-{\frac {1}{2}}S_{5}-{\frac {1}{2}}S_{6}-{\frac {1}{4}}S_{7},}
( 2 ) ln ⁡ ⁡ 2 = 2 R e ( I ) = S 1 − − 1 2 S 3 − − 1 2 S 4 − − 1 4 S 5 + 1 8 S 7 + 1 8 S 8 . {\displaystyle (2)\qquad \ln 2=2\ \mathrm {Re} (I)=S_{1}-{\frac {1}{2}}S_{3}-{\frac {1}{2}}S_{4}-{\frac {1}{4}}S_{5}+{\frac {1}{8}}S_{7}+{\frac {1}{8}}S_{8}.}

Mais ln 2 s'exprime par ailleurs plus directement en fonction des Sn :

( 3 ) ln ⁡ ⁡ 2 = [ − − ln ⁡ ⁡ ( 2 − − y 2 ) ] 0 1 = ∫ ∫ 0 1 y 1 − − y 2 / 2 d y = ∑ ∑ k ≥ ≥ 0 1 2 ( k + 1 ) 2 k = S 2 + 1 2 S 4 + 1 4 S 6 + 1 8 S 8 . {\displaystyle {\begin{aligned}(3)\quad \ln 2&=[-\ln(2-y^{2})]_{0}^{1}=\int _{0}^{1}{\frac {y}{1-y^{2}/2}}\mathrm {d} y=\sum _{k\geq 0}{\frac {1}{2(k+1)2^{k}}}\\&=S_{2}+{\frac {1}{2}}S_{4}+{\frac {1}{4}}S_{6}+{\frac {1}{8}}S_{8}.\end{aligned}}}

Les expressions (2) et (3) donnent donc, par soustraction des membres de droite, une relation entre les Sn :

( 4 ) 0 = 2 S 1 − − 2 S 2 − − S 3 − − 2 S 4 − − 1 2 S 5 − − 1 2 S 6 + 1 4 S 7 . {\displaystyle (4)\qquad 0=2S_{1}-2S_{2}-S_{3}-2S_{4}-{\frac {1}{2}}S_{5}-{\frac {1}{2}}S_{6}+{\frac {1}{4}}S_{7}.}

En multipliant le membre de droite de (4) par 1 + 4r et en ajoutant ce produit au membre de droite de (1), on obtient l'égalité (0) annoncée.

Exploitation de la formule pour calculer les chiffres après la virgule de π

Le but est de calculer le N-ième chiffre après la virgule de π en base 16.

Déjà, on remarque que le (N + 1)-ième chiffre après la virgule de π en base 16 est le même que le 1er chiffre après la virgule de 16Nπ. En effet, comme en base 10, multiplier un nombre en base 16 par 16 permet de décaler la virgule d'un rang vers la droite. En multipliant un nombre par 16N, la virgule est donc décalée de N rangs vers la droite. Ainsi, il suffit de calculer le premier chiffre de 16N π, égal par la formule BBP à :

16

N

π

π

=



k

=

0



16

N



k

(

4

8

k

+

1



2

8

k

+

4



1

8

k

+

5



1

8

k

+

6

)

.

{\displaystyle 16^{N}\pi =\sum _{k=0}^{\infty }16^{N-k}\left({\frac {4}{8k+1}}-{\frac {2}{8k+4}}-{\frac {1}{8k+5}}-{\frac {1}{8k+6}}\right).}

Mais calculer les premiers chiffres derrière la virgule de ce nombre n'est pas si simple, pour deux raisons :

• d'abord, ce nombre étant très grand ( N ≫ ≫ 1 ) , {\textstyle (N\gg 1),} cela demande d'effectuer des calculs sur des nombres très grands ;
• ensuite, parce que cette sommation est infinie.

On pose S N ( a ) = ∑ ∑ k = 0 ∞ ∞ 16 N − − k 8 k + a . {\displaystyle S_{N}(a)=\sum _{k=0}^{\infty }{\frac {16^{N-k}}{8k+a}}.}

Le calcul des premiers chiffres de SN(a) permettra d'obtenir ceux de 16N π ( a = 1 , 4..6 ) {\textstyle (a=1,4..6)} , par la relation :

16

N

π

π

=

4

S

N

(

1

)



2

S

N

(

4

)



S

N

(

5

)



S

N

(

6

)

.

{\displaystyle 16^{N}\pi =4\ S_{N}(1)-2\ S_{N}(4)-S_{N}(5)-S_{N}(6).}

On découpe la somme SN(a) en deux :

S

N

(

a

)

=



k

=

0

N



1

16

N



k

8

k

+

a

+



k

=

N



16

N



k

8

k

+

a

=

A

N

(

a

)

+

B

N

(

a

)

{\displaystyle S_{N}(a)=\sum _{k=0}^{N-1}{\frac {16^{N-k}}{8k+a}}+\sum _{k=N}^{\infty }{\frac {16^{N-k}}{8k+a}}=A_{N}(a)+B_{N}(a)}

et on calcule AN(a) et BN(a) indépendamment.

Calcul de B N ( a )

B N ( a ) = ∑ ∑ k = N ∞ ∞ 16 N − − k 8 k + a {\displaystyle B_{N}(a)=\sum _{k=N}^{\infty }{\frac {16^{N-k}}{8k+a}}}

Bien que ce soit une somme infinie, ce terme est très simple à calculer, car on remarque que ses termes b N ( k , a ) {\textstyle b_{N}(k,a)} deviennent vite très petits et on ne cherche que les premiers chiffres.

• En effet, le premier terme de la somme est : b N = 1 8 N + a . {\textstyle b_{N}={\frac {1}{8N+a}}.} Comme on cherche le N-ième chiffre derrière la virgule de π (N = 1 000 000 000 par exemple), le premier terme bN est très inférieur à 1.
• De plus, chaque terme suivant a un zéro de plus derrière la virgule que le précédent, car pour kN, bk > 16 bk+1 :

b

k

b

k

+

1

=

16

N



k

16

N



(

k

+

1

)

8

(

k

+

1

)

+

a

8

k

+

a

=

16

(

1

+

8

8

k

+

a

)



16

+

.

{\displaystyle {\frac {b_{k}}{b_{k+1}}}={\frac {16^{N-k}}{16^{N-(k+1)}}}{\frac {8(k+1)+a}{8k+a}}=16\left(1+{\frac {8}{8k+a}}\right)\longrightarrow 16^{+}.}

Finalement, la somme BN(a) est de la forme (au pire) :

B N = 0 , ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ . . . {\displaystyle B_{N}\ =\ \ 0,**********\ .\ .\ .}

+ 0 , 0 ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ . . . {\displaystyle +\ 0,0*********\ .\ .\ .}

+ 0 , 00 ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ . . . {\displaystyle +\ 0,00********\ .\ .\ .}

+ 0 , 000 ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ . . . {\displaystyle +\ 0,000********\ .\ .\ .}

Donc pour obtenir BN(a) avec une précision de P chiffres derrière la virgule, il suffit de calculer les P premiers termes de la somme, plus les quelques suivants pour éviter les problèmes de retenues qui peuvent éventuellement apparaître.

Il suffit donc de calculer : B N ′ ( a ) = ∑ ∑ k = N N + P + 10 16 N − − k 8 k + a = 16 0 8 N + a + 16 − − 1 8 ( N + 1 ) + a + . . . + 16 − − ( P + 10 ) 8 ( N + P + 10 ) + a . {\displaystyle B_{N}'(a)=\sum _{k=N}^{N+P+10}{\frac {16^{N-k}}{8k+a}}={\frac {16^{0}}{8N+a}}+{\frac {16^{-1}}{8(N+1)+a}}+...+{\frac {16^{-(P+10)}}{8(N+P+10)+a}}.}

Cette somme n'étant composée que d'un petit nombre de termes P + 11 {\textstyle P+11} (nombre constant), son temps de calcul est négligeable pour un ordinateur.

Calcul de A N ( a )

A N ( a ) = ∑ ∑ k = 0 N − − 1 16 N − − k 8 k + a = 16 N a + 16 N − − 1 8 + a + . . . + 16 1 8 ( N − − 1 ) + a . {\displaystyle A_{N}(a)=\sum _{k=0}^{N-1}{\frac {16^{N-k}}{8k+a}}={\frac {16^{N}}{a}}+{\frac {16^{N-1}}{8+a}}+...+{\frac {16^{1}}{8(N-1)+a}}.}

Le problème pour calculer AN(a) est que les premiers termes sont extrêmement grands (N chiffres en base 16 devant la virgule). Néanmoins, comme on ne cherche que les premiers chiffres derrière la virgule, peu importe la partie entière, aussi grande qu'elle soit. On peut donc s'en « débarrasser » en utilisant l'arithmétique modulaire.

Toute la difficulté se réduit donc à trouver la partie fractionnelle de 16 N − − k 8 k + a . {\displaystyle {\frac {16^{N-k}}{8k+a}}.}

Pour cela, on effectue la division euclidienne de 16N-k par 8k+a :

∃ ∃ q ∈ ∈ Z , ∃ ∃ r < 8 k + a , 16 N − − k = q ( 8 k + a ) + r {\displaystyle \exists q\in \mathbb {Z} ,\exists r<8k+a,\ 16^{N-k}=q(8k+a)+r}

Donc 16 N − − k 8 k + a = q + r 8 k + a {\displaystyle {\frac {16^{N-k}}{8k+a}}=q+{\frac {r}{8k+a}}}

r 8 k + a {\textstyle {\frac {r}{8k+a}}} est inférieur à 1, donc c'est la partie fractionnelle de 16 N − − k 8 k + a . {\textstyle {\frac {16^{N-k}}{8k+a}}.}

Et r 8 k + a = 16 N − − k mod ( 8 k + a ) 8 k + a . {\textstyle {\frac {r}{8k+a}}={\frac {16^{N-k}{\bmod {(}}8k+a)}{8k+a}}.}

Il suffit donc de calculer :

A N ′ ( a ) = ∑ ∑ k = 0 N − − 1 16 N − − k mod ( 8 k + a ) 8 k + a = 16 N mod ( a ) a + 16 N − − 1 mod ( 8 + a ) 8 + a + 16 N − − 2 mod ( 16 + a ) 16 + a + . . . {\displaystyle A_{N}'(a)=\sum _{k=0}^{N-1}{\frac {16^{N-k}{\bmod {(}}8k+a)}{8k+a}}={\frac {16^{N}{\bmod {(}}a)}{a}}+{\frac {16^{N-1}{\bmod {(}}8+a)}{8+a}}+{\frac {16^{N-2}{\bmod {(}}16+a)}{16+a}}+...}

En utilisant la méthode d'exponentiation rapide, 16N-k mod (8k+a) se calcule rapidement (temps d'exécution en O(log2(N-k)) ). Temps totale pour N ≫ ≫ 1 {\textstyle N\gg 1} termes est O ( log 2 ⁡ ⁡ N ) + O ( log 2 ⁡ ⁡ ( N − − 1 ) ) + O ( log 2 ⁡ ⁡ ( N − − 2 ) ) + ⋯ ⋯ {\textstyle O(\log _{2}N)+O(\log _{2}(N-1))+O(\log _{2}(N-2))+\cdots }

Conclusion

Finalement, pour obtenir les premiers chiffres de π en base 16 (ou 2), il faut calculer les premiers chiffres de :

π π N = 4 S N ′ ( 1 ) − − 2 S N ′ ( 4 ) − − S N ′ ( 5 ) − − S N ′ ( 6 ) {\displaystyle \pi _{N}=4\ S_{N}'(1)-2\ S_{N}'(4)-S_{N}'(5)-S_{N}'(6)}

avec S N ′ ( a ) = ∑ ∑ k = 0 N − − 1 16 N − − k mod ( 8 k + a ) 8 k + a + ∑ ∑ k = N N + P + 10 16 N − − k 8 k + a . {\displaystyle S_{N}'(a)=\sum _{k=0}^{N-1}{\frac {16^{N-k}{\bmod {(}}8k+a)}{8k+a}}+\sum _{k=N}^{N+P+10}{\frac {16^{N-k}}{8k+a}}.}

Complexité de cette méthode

Pour calculer le n-ième chiffre après la virgule de π en base 16 (et donc le 4n-ième chiffre en base 2) :

Complexité temporelle

BN'(a) se calcule en temps constant (O(1)). La complexité du calcul de SN' est donc la même que la complexité du calcul de AN'(a).
AN'(a) : en utilisant la méthode d'exponentiation rapide, ses termes se calculent en O(log2(n)) multiplications sur des entiers de taille log2(n). En notant M(k) la complexité de la multiplication de deux entiers de taille k, la complexité est donc O(log2(n)M(log2(n)) ). Finalement, la somme des n termes, AN'(a), se calcule en temps O(n log2(n)M(log2(n)) ). Même en utilisant l'algorithme de multiplication naïf plutôt que l'algorithme de Karatsuba ou la transformée de Fourier rapide, on obtient une complexité quasi linéaire de O(n log2(n)3).

Complexité spatiale

Le calcul de BN'(a) s'effectue en espace constant (somme d'un nombre fixé de termes, avec un nombre fixé de chiffres significatifs). Le calcul de AN'(a) nécessite d'effectuer des calculs modulo 8k+a, c'est-à-dire de manipuler des nombres de taille log2(k) avec k < N. À chaque étape de l'algorithme, on manipule un nombre constant de tels nombres : la complexité en espace du calcul de AN'(a) est donc O(log2(n)). L'algorithme total utilise donc un espace logarithmique.

Formules dérivées

Simon Plouffe (1995)

Formule originale : π π = ∑ ∑ k = 0 ∞ ∞ 1 16 k ( 4 8 k + 1 − − 2 8 k + 4 − − 1 8 k + 5 − − 1 8 k + 6 ) {\displaystyle \pi \ =\ \sum _{k=0}^{\infty }{\frac {1}{16^{k}}}\left({\frac {4}{8k+1}}-{\frac {2}{8k+4}}-{\frac {1}{8k+5}}-{\frac {1}{8k+6}}\right)}

∀ ∀ r ∈ ∈ C π π = ∑ ∑ k = 0 ∞ ∞ 1 16 k ( 4 + 8 r 8 k + 1 − − 8 r 8 k + 2 − − 4 r 8 k + 3 − − 2 + 8 r 8 k + 4 − − 1 + 2 r 8 k + 5 − − 1 + 2 r 8 k + 6 + r 8 k + 7 ) {\displaystyle \forall r\in \mathbb {C} \ \ \ \ \pi \ =\ \sum _{k=0}^{\infty }{\frac {1}{16^{k}}}\left({\frac {4+8r}{8k+1}}-{\frac {8r}{8k+2}}-{\frac {4r}{8k+3}}-{\frac {2+8r}{8k+4}}-{\frac {1+2r}{8k+5}}-{\frac {1+2r}{8k+6}}+{\frac {r}{8k+7}}\right)}

π π 2 = ∑ ∑ k = 0 ∞ ∞ ( − − 1 ) k 8 k ( 4 6 k + 1 + 1 6 k + 2 + 1 6 k + 3 ) {\displaystyle \pi {\sqrt {2}}\ =\ \sum _{k=0}^{\infty }{\frac {(-1)^{k}}{8^{k}}}\left({\frac {4}{6k+1}}+{\frac {1}{6k+2}}+{\frac {1}{6k+3}}\right)}

π π 2 = ∑ ∑ k = 0 ∞ ∞ 1 16 k ( 16 ( 8 k + 1 ) 2 − − 16 ( 8 k + 2 ) 2 − − 8 ( 8 k + 3 ) 2 − − 16 ( 8 k + 4 ) 2 − − 4 ( 8 k + 5 ) 2 − − 4 ( 8 k + 6 ) 2 − − 2 ( 8 k + 7 ) ) 2 ) {\displaystyle \pi ^{2}\ =\ \sum _{k=0}^{\infty }{\frac {1}{16^{k}}}\left({\frac {16}{(8k+1)^{2}}}-{\frac {16}{(8k+2)^{2}}}-{\frac {8}{(8k+3)^{2}}}-{\frac {16}{(8k+4)^{2}}}-{\frac {4}{(8k+5)^{2}}}-{\frac {4}{(8k+6)^{2}}}-{\frac {2}{(8k+7))^{2}}}\right)}

π π 2 = 9 8 ∑ ∑ k = 0 ∞ ∞ 1 64 k ( 16 ( 6 k + 1 ) 2 − − 24 ( 6 k + 2 ) 2 − − 8 ( 6 k + 3 ) 2 − − 6 ( 6 k + 4 ) 2 − − 1 ( 6 k + 5 ) 2 ) {\displaystyle \pi ^{2}\ =\ {\frac {9}{8}}\ \sum _{k=0}^{\infty }{\frac {1}{64^{k}}}\left({\frac {16}{(6k+1)^{2}}}-{\frac {24}{(6k+2)^{2}}}-{\frac {8}{(6k+3)^{2}}}-{\frac {6}{(6k+4)^{2}}}-{\frac {1}{(6k+5)^{2}}}\right)}

π π 2 = 2 27 ∑ ∑ k = 0 ∞ ∞ 1 729 k ( 243 ( 12 k + 1 ) 2 − − 405 ( 12 k + 2 ) 2 − − 81 ( 12 k + 4 ) 2 − − 27 ( 12 k + 5 ) 2 − − 72 ( 12 k + 6 ) 2 − − 9 ( 12 k + 7 ) 2 − − 9 ( 12 k + 8 ) 2 − − 5 ( 12 k + 10 ) 2 + 1 ( 12 k + 11 ) 2 ) {\displaystyle \pi ^{2}\ =\ {\frac {2}{27}}\ \sum _{k=0}^{\infty }{\frac {1}{729^{k}}}\left({\frac {243}{(12k+1)^{2}}}-{\frac {405}{(12k+2)^{2}}}-{\frac {81}{(12k+4)^{2}}}-{\frac {27}{(12k+5)^{2}}}-{\frac {72}{(12k+6)^{2}}}-{\frac {9}{(12k+7)^{2}}}-{\frac {9}{(12k+8)^{2}}}-{\frac {5}{(12k+10)^{2}}}+{\frac {1}{(12k+11)^{2}}}\right)}

Adamchick et Wagon (1997)

π π = ∑ ∑ i = 0 ∞ ∞ ( − − 1 ) i 4 i ( 2 4 i + 1 + 2 4 i + 2 + 1 4 i + 3 ) {\displaystyle \pi \ =\ \sum _{i=0}^{\infty }{\frac {(-1)^{i}}{4^{i}}}\left({\frac {2}{4i+1}}+{\frac {2}{4i+2}}+{\frac {1}{4i+3}}\right)}

Fabrice Bellard

π π = 1 64 ∑ ∑ i = 0 ∞ ∞ ( − − 1 ) i 2 10 i ( − − 32 4 i + 1 − − 1 4 i + 3 + 256 10 i + 1 − − 64 10 i + 3 − − 4 10 i + 5 − − 4 10 i + 7 + 1 10 i + 9 ) {\displaystyle \pi \ =\ {\frac {1}{64}}\ \sum _{i=0}^{\infty }{\frac {(-1)^{i}}{2^{10i}}}\left(-{\frac {32}{4i+1}}-{\frac {1}{4i+3}}+{\frac {256}{10i+1}}-{\frac {64}{10i+3}}-{\frac {4}{10i+5}}-{\frac {4}{10i+7}}+{\frac {1}{10i+9}}\right)}

Géry Huvent (2001)

π π = 1 128 ∑ ∑ i = 0 ∞ ∞ 1 2 12 i ( 768 24 i + 3 + 512 24 i + 4 + 128 24 i + 6 − − 16 24 i + 12 − − 16 24 i + 14 − − 12 24 i + 15 + 2 24 i + 20 − − 1 24 i + 22 ) {\displaystyle \pi \ =\ {\frac {1}{128}}\ \sum _{i=0}^{\infty }{\frac {1}{2^{12i}}}\left({\frac {768}{24i+3}}+{\frac {512}{24i+4}}+{\frac {128}{24i+6}}-{\frac {16}{24i+12}}-{\frac {16}{24i+14}}-{\frac {12}{24i+15}}+{\frac {2}{24i+20}}-{\frac {1}{24i+22}}\right)}
π π = 1 256 ∑ ∑ i = 0 ∞ ∞ 1 2 12 i ( 2048 24 i + 2 − − 1024 24 i + 4 + 192 24 i + 9 + 128 24 i + 10 + 32 24 i + 12 − − 4 24 i + 18 − − 4 24 i + 20 − − 3 24 i + 21 ) {\displaystyle \pi \ =\ {\frac {1}{256}}\ \sum _{i=0}^{\infty }{\frac {1}{2^{12i}}}\left({\frac {2048}{24i+2}}-{\frac {1024}{24i+4}}+{\frac {192}{24i+9}}+{\frac {128}{24i+10}}+{\frac {32}{24i+12}}-{\frac {4}{24i+18}}-{\frac {4}{24i+20}}-{\frac {3}{24i+21}}\right)}
π π 3 = 1 16 ∑ ∑ i = 0 ∞ ∞ ( − − 1 ) i 2 10 i ( 32 ( 4 i + 1 ) 3 + 8 ( 4 i + 2 ) 3 + 1 ( 4 i + 3 ) 3 ) {\displaystyle \pi ^{3}\ =\ {\frac {1}{16}}\ \sum _{i=0}^{\infty }{\frac {(-1)^{i}}{2^{10i}}}\left({\frac {32}{(4i+1)^{3}}}+{\frac {8}{(4i+2)^{3}}}+{\frac {1}{(4i+3)^{3}}}\right)} + 5 2 ∑ ∑ i = 0 ∞ ∞ ( − − 1 ) i 2 6 i ( 32 ( 12 i + 1 ) 3 − − 192 ( 12 i + 2 ) 3 + 88 ( 12 i + 3 ) 3 − − 8 ( 12 i + 5 ) 3 + 84 ( 12 i + 6 ) 3 − − 4 ( 12 i + 7 ) 3 + 11 ( 12 i + 9 ) 3 − − 12 ( 12 i + 10 ) 3 + 1 ( 12 i + 11 ) 3 ) {\displaystyle \ +\ {\frac {5}{2}}\ \sum _{i=0}^{\infty }{\frac {(-1)^{i}}{2^{6i}}}\left({\frac {32}{(12i+1)^{3}}}-{\frac {192}{(12i+2)^{3}}}+{\frac {88}{(12i+3)^{3}}}-{\frac {8}{(12i+5)^{3}}}+{\frac {84}{(12i+6)^{3}}}-{\frac {4}{(12i+7)^{3}}}+{\frac {11}{(12i+9)^{3}}}-{\frac {12}{(12i+10)^{3}}}+{\frac {1}{(12i+11)^{3}}}\right)}

π π 4 = 27 164 ∑ ∑ i = 0 ∞ ∞ 1 2 12 i ( 2048 ( 24 i + 1 ) 4 − − 38912 ( 24 i + 2 ) 4 + 81920 ( 24 i + 3 ) 4 − − 2048 ( 24 i + 4 ) 4 − − 512 ( 24 i + 5 ) 4 − − 23552 ( 24 i + 6 ) 4 + 256 ( 24 i + 7 ) 4 − − 27648 ( 24 i + 8 ) 4 − − 10240 ( 24 i + 9 ) 4 {\displaystyle \pi ^{4}\ =\ {\frac {27}{164}}\ \sum _{i=0}^{\infty }{\frac {1}{2^{12i}}}\left({\frac {2048}{(24i+1)^{4}}}-{\frac {38912}{(24i+2)^{4}}}+{\frac {81920}{(24i+3)^{4}}}-{\frac {2048}{(24i+4)^{4}}}-{\frac {512}{(24i+5)^{4}}}-{\frac {23552}{(24i+6)^{4}}}+{\frac {256}{(24i+7)^{4}}}-{\frac {27648}{(24i+8)^{4}}}-{\frac {10240}{(24i+9)^{4}}}\right.} − − 2432 ( 24 i + 10 ) 4 − − 64 ( 24 i + 11 ) 4 − − 3584 ( 24 i + 12 ) 4 − − 32 ( 24 i + 13 ) 4 − − 608 ( 24 i + 14 ) 4 − − 1280 ( 24 i + 15 ) 4 − − 1728 ( 24 i + 16 ) 4 + 8 ( 24 i + 17 ) 4 − − 368 ( 24 i + 18 ) 4 − − 4 ( 24 i + 19 ) 4 {\displaystyle -{\frac {2432}{(24i+10)^{4}}}-{\frac {64}{(24i+11)^{4}}}-{\frac {3584}{(24i+12)^{4}}}-{\frac {32}{(24i+13)^{4}}}-{\frac {608}{(24i+14)^{4}}}-{\frac {1280}{(24i+15)^{4}}}-{\frac {1728}{(24i+16)^{4}}}+{\frac {8}{(24i+17)^{4}}}-{\frac {368}{(24i+18)^{4}}}-{\frac {4}{(24i+19)^{4}}}} − − 8 ( 24 i + 20 ) 4 + 160 ( 24 i + 21 ) 4 − − 38 ( 24 i + 22 ) 4 + 1 ( 24 i + 23 ) 4 ) {\displaystyle -\left.{\frac {8}{(24i+20)^{4}}}+{\frac {160}{(24i+21)^{4}}}-{\frac {38}{(24i+22)^{4}}}+{\frac {1}{(24i+23)^{4}}}\right)} cite-ref-2[2]

Les records

Pour comparaison, le record de calcul de toutes les décimales de π est, en 2016, de 22 600 milliards de décimales (soit environ 70 000 milliards de chiffres binaires).

• 7 octobre 1996 (Fabrice Bellard) : 400 milliardième chiffre en base 2
• septembre 1997 (Fabrice Bellard) : 1 000 milliardième chiffre en base 2
• février 1999 (Colin Percival) : 40 000 milliardième chiffre en base 2
2001 : 4 000 000 milliardième chiffre en base 2

Calcul en base 10

Actuellement, aucune formule réellement efficace n'a été découverte pour calculer le n-ième chiffre de π en base 10. Simon Plouffe a mis au point en décembre 1996, à partir d'une très ancienne série de calcul de π basée sur les coefficients du binôme de Newton, une méthode pour calculer les chiffres en base 10, mais sa complexité en O(n3 log2(n)) la rendait en pratique inutilisable. Fabrice Bellard a bien amélioré l'algorithme pour atteindre une complexité en O(n2), mais cela n'est pas suffisant pour concurrencer les méthodes classiques de calcul de toutes les décimales.

Formule de type BBP

Par extension des formules BBP qui permettent un calcul de π, Barsky, V. Muñoz et R. Pérez-Marco ont généralisé la forme des séries à la base des formules BBP en posantcite-ref-3[3]cite-ref-4[4]:

B B P ( d , n , b , A ) = ∑ ∑ p = 0 + ∞ ∞ 1 b p ∑ ∑ j = 1 n a j ( k n + j ) p {\displaystyle \mathrm {BBP} (d,n,b,A)=\sum _{p=0}^{+\infty }{\frac {1}{b^{p}}}\sum _{j=1}^{n}{\frac {a_{j}}{(kn+j)^{p}}}}

où on désigne par d le degré, b la base, n le nombre et A, un vecteur de longueur n bien choisi.

Ainsi, la formule originelle de Plouffe s'obtient avec d = 1, b = 16, n = 8, A = (4,0,0,−2,−1,−1,0,0).

Le principe de la généralisation repose sur l'utilisation du développement en série entière de la fonction z → log(1 + z) appliqué à des valeurs complexes bien choisies, ce qui permet de calculer des valeurs de puissances de π, mais aussi des valeurs des fonctions zêta de Riemann, bêta et êta de Dirichlet ou dilogarithmecite-ref-5[5].

Références

cite-note-11. baileyborweinplouffe1997david-h-baileypeter-b-borweinsimon-plouffe1997(en) David H. Bailey, Peter B. Borwein et Simon Plouffe, « On the Rapid Computation of Various Poly-logarithmic Constants », Math. Comp., vol. 66, no 218,‎ 1997, p. 903-13 (DOI 10.1090/S0025-5718-97-00856-9, MR 1415794).
cite-note-22. huvent2001g-ry-huvent2001Géry Huvent, « Formules BBP », 2001.
cite-note-33. (en) Daniel Barsky, Vicente Muñoz et Ricardo Pérez-Marco, « On the genesis of BBP formulas », 23 juin 2019.
cite-note-44. kristensenmathiasen2023simon-kristensenoskar-mathiasen2023(en) Simon Kristensen et Oskar Mathiasen, « BBP-type formulas — An elementary approach », Journal of Number Theory, vol. 244,‎ mars 2023, p. 251-263 (DOI 10.1016/j.jnt.2022.09.001)
cite-note-55. (en) D. J. Broadhurst, « Polylogarithmic ladders, hypergeometric series and the ten millionth digits of ζ(3) and ζ(5) », .

Bibliographie

• arndthaenel2006j-rg-arndtchristoph-haenel2006Jörg Arndt et Christoph Haenel (trad. de l'allemand par Henri Lemberg et François Guénard), À la Poursuite de π [« Pi. Algorithmen, Computer, Arithmetik »], Paris, Vuibert,‎ mars 2006 (1re éd. 1998), 273 p. (ISBN 978-2-7117-7170-7), chap. 10 (« L'algorithme BBP »).
• adamchikwagon1997victor-adamchikstan-wagon1997(en) Victor Adamchik et Stan Wagon, « A Simple Formula for π », The American Mathematical Monthly, vol. 104, no 9,‎ 1997, p. 852–55 (DOI 10.1080/00029890.1997.11990729)

Liens externes

weissteineric-w-weisstein(en) Eric W. Weisstein, « BBP Formula », sur MathWorld

• Portail des mathématiques